%% Directory

% Matlab needs to be configured such that it can find the fun_***.m files
% one way to do this is 'addpath'
addpath 'LOCATION OF THE MATLAB .m FILE FOLDER HERE'


%% Exogenous parameters
clear all

T=85;
t=58;
t_retstart = 62; 
t_ret = 66;
d=2; 
beta=0.98;
R_base = 1/beta;
R_shocked = R_base - 0.001;     % 10 BASIS POINT SHOCK

probgiving = 0.12;  


logshifter = 10/1000000;

grossinc = 0.77;  % everything in millions

c_hat = 0.93*grossinc*(1-0.40);


netinc = grossinc*(1-0.40); % inc tax
initialwealth = 4.23;

s_vector = (1:1:T)';

p = NaN(T,1);
p_base = NaN(T,1);
p_base(t:t+5) = [0.72;0.73;0.73;0.75;0.76;0.77];
p_base(t+6:T) = 0.77;

w=NaN(T,1);
w(t) = netinc + initialwealth;
w(t+1:t_retstart) = netinc;
aidur = 0.5*1/(t_ret-t_retstart); % annual income decrease until retirement
w(t_retstart:T) = netinc - min(netinc*aidur*( s_vector(t_retstart:T)-t_retstart), netinc*0.5);

g_hat = 5225/1000000;

% to be filled in
g = NaN(T,1);
c = NaN(T,1);

p = p_base;

varepsilon = 0.44;

empirical_cond_TE = -0.2603*0.1;

%% Using conditional giving response to inform EIS

sigmavector = [0.03;0.04;0.05;0.055;0.058;0.059;0.06;0.061;0.062;0.063;0.06303;0.06305;0.064;0.065;0.07;0.08;0.09;(0.1:0.01:0.41)']; % numerical issues if pulled further
g_cond_control = NaN(size(sigmavector));
g_cond_treat   = NaN(size(sigmavector));

for isigma=1:size(sigmavector,1)
    kappac = c_hat^(-1/sigmavector(isigma))*g_hat^(1/varepsilon)*p(t);
    disp(num2str(isigma))
    % Control group, "1" arg gives us giving amount as output
    g_cond_control(isigma) = fun_optimize_wEC(1, kappac,sigmavector(isigma),varepsilon,beta, R_base, t, T, 0, p, s_vector, w);
        
    % Treatment group
    g_cond_treat(isigma)   = fun_optimize_wEC(1, kappac,sigmavector(isigma),varepsilon,beta, R_shocked, t, T, 0, p, s_vector, w);
end

cond_treateffect = log(logshifter+g_cond_treat)-log(logshifter+g_cond_control);

%% Finding the EIS
fun_kappac = @(x) c_hat^(-1/ x )*g_hat^(1/varepsilon)*p(t)
fun_TE = @(x)    double( log(logshifter+fun_optimize_wEC(1, fun_kappac(x),  x  ,varepsilon,beta, R_shocked, t, T, 0, p, s_vector, w)) ...
                       - log(logshifter+fun_optimize_wEC(1, fun_kappac(x),  x  ,varepsilon,beta, R_base,    t, T, 0, p, s_vector, w)));

fun_TE_diff_sq = @(x) double( (fun_TE(x)- empirical_cond_TE)^2 );

options = optimoptions('fmincon','GradObj','off','OptimalityTolerance', 1e-10,'StepTolerance',1e-10,'Display','iter-detailed');
    
sigma_est  = fmincon(fun_TE_diff_sq,0.007,-1,-0.02,[],[],[],[],[],options )


%% Giving Response for Different Values of the MEC
Nh = 4999;

sigma = sigma_est; 
kappac = c_hat^(-1/sigma)*g_hat^(1/varepsilon)*p(t);

MECvector = [10000;20000;22500;25000;27500;30000;31000;32000;33000;35000;40000;45000;50000;55000;60000;65000;70000;75000;80000;85000;90000;95000;100000]/1000000;
dP = NaN(size(MECvector));
LNsigmas = NaN(size(MECvector));
kappa_stars = NaN(size(MECvector));
kappavector = NaN(size(MECvector,1),Nh);
psivector1 = NaN(size(MECvector,1),1);
psivector2 = NaN(size(MECvector,1),1);
utildiff = NaN(size(MECvector,1),1);
P_control = NaN(size(MECvector,1),1);
for iMEC = 1:size(MECvector,1)

    % Find kappa such that ppl are indifferent between giving or not
    % because in order to match empirical P(g>0), need to know what
    % the kappa distribution should look like
    
    % @ function that calls the fun_to_optimize_wEC function
    % --> first arg=2 so that function returns U(giving>0) - U(giving = 0, 
    %     with MEC added to nongiver's endowment)
    fun_to_minimize = @(x) fun_optimize_wEC(2,abs(x),sigma,varepsilon,beta, R_base, t, T, MECvector(iMEC), p, s_vector, w)^2;
    
    % Minimize the square of the utility difference
    options = optimoptions('fmincon','GradObj','off');
 
    [kappa_stars(iMEC),utildiff(iMEC)] = fmincon(fun_to_minimize,kappac,[],[],[],[],[],[],[],options ); % init guess = kappac
    kappa_stars(iMEC)=abs(kappa_stars(iMEC));

    % Find distribution parameters
    options = optimoptions('fmincon','GradObj','off','OptimalityTolerance', 1e-10,'StepTolerance',1e-20);
    
    A = [0 0; 0 -1];
    kappa_star_ = kappa_stars(iMEC);
    sol  = fmincon(@(x) fun_to_minimize_to_fit_lognormal(kappac,kappa_star_,probgiving,x),[log(kappac);1],A,[0;0],[],[],[],[],[],options )';
    psivector1(iMEC) = sol(1);
    psivector2(iMEC) = sol(2);
    
    % Define the resulting probability distribution as an object
    pd = makedist('LogNormal','mu',psivector1(iMEC),'sigma',psivector2(iMEC));
    
    LNsigmas(iMEC) = psivector2(iMEC);
    
    g_control = NaN(Nh,1);
    g_treat   = NaN(Nh,1);

    parfor(h = 1:Nh,25)
        disp(['========> iMEC = ' num2str(iMEC) ' / ' num2str(size(MECvector,1)) '   .....  h = ' num2str(h) ' / ' num2str(Nh)])
        kappavector(iMEC,h)=icdf('LogNormal',h/(Nh+1), psivector1(iMEC),psivector2(iMEC)  );

        % Control group, "1" arg gives us giving amount as output
        g1(h) = fun_optimize_wEC(1, kappavector(iMEC,h),sigma,varepsilon,beta, R_base, t, T, MECvector(iMEC), p, s_vector, w);
    
        % Treatment group
        g2(h) = fun_optimize_wEC(1, kappavector(iMEC,h),sigma,varepsilon,beta, R_shocked, t, T, MECvector(iMEC), p, s_vector, w);
    end
    P_control(iMEC) = mean(g(:,1)>0);
    dP(iMEC) = mean((g2(:)>0) - (g1(:)>0));

end


%% We also want to have the amount of giving of the marginal untreated giver
% get this by using the marginal type (kappastar) and give them slightly lower 
% lower MEC.
marg_pos_amount = NaN(size(MECvector));
lifetimegiving_marggiver = NaN(size(MECvector));
for iMEC = 1:size(MECvector,1)
    marg_pos_amount(iMEC) = fun_optimize_wEC(1, kappa_stars(iMEC),sigma,varepsilon,beta, R_base, t, T, MECvector(iMEC)*(1-0.0001), p, s_vector, w);
    
    % also what is the PV of life-time giving?
    gvec = fun_optimize_wEC(3, kappa_stars(iMEC),sigma,varepsilon,beta, R_base, t, T, MECvector(iMEC)*(1-0.0001), p, s_vector, w) ;
    lifetimegiving_marggiver(iMEC) = ones(T-t+1,1)'*R_base^(-1)*gvec(t:T);
end

MEC_divby_lifetimegiving_marggiver = MECvector./lifetimegiving_marggiver

% our solution is at iMEC=10. What is MEC divided by Present Value of
% Life-time giving?
MEC_divby_lifetimegiving_marggiver(10)


%% print results for copying into stata for plotting


[MECvector dP]

% Paste this in below '	0.000	   .' and before 'end'
% should take up lines 9 through 32 in wdon_stata_plot_extensivemargin.do

%% Counterfactual: How much would giving increase if MEC=0
% These numbers are referenced in text
iMEC=10;
g1_noEC=NaN(Nh,1);

    parfor(h = 1:Nh,25)
        disp(['========> iMEC = ' num2str(iMEC) ' / ' num2str(size(MECvector,1)) '   .....  h = ' num2str(h) ' / ' num2str(Nh)])
       
        % Control group, "1" arg gives us giving amount as output
        g1_noEC(h) = fun_optimize_wEC(1, kappavector(iMEC,h),sigma,varepsilon,beta, R_base, t, T, 0, p, s_vector, w);
    end

% mean amount of giving among nongivers 
[kappavector(iMEC,:) < kappa_stars(iMEC)]*g1_noEC/sum([kappavector(iMEC,:) < kappa_stars(iMEC)])

% Effect on mean, divided by existing unconditional mean
[kappavector(iMEC,:) < kappa_stars(iMEC)]*g1_noEC/sum([kappavector(iMEC,:) < kappa_stars(iMEC)])/ mean( g1_noEC' .* [kappavector(iMEC,:) >= kappa_stars(iMEC)] )





